Example 1: Single Data Analysis Unit

Author: Brian D. Gerber

Date: 2023-05-10

The primary objective of the Diel.Niche package is to be able to use observations of wild animals to estimate support for a set of diel niche hypotheses. The fundamental data analysis unit is a set of three frequencies (y): the number of observations of a species during twilight, daytime, and nighttime (in that order). In this vignette, we will demonstrate the basic functionality of the package. This includes 1) model comparison, 2) parameter estimation, 3) posterior predictive check, and 4) results visualization.

We will use data available in the package (‘diel.data’) provided by the Urban Wildlife Information Network (https://www.urbanwildlifeinfo.org/). The specific data are camera trap detections of the Virginia Opossum (Didelphis virginiana) during the winter of 2018. The available data are aggregated independent counts from 131 camera locations in Chicago, Illinois USA. Our objective is to evaluate the support for the traditional diel hypotheses of the Virginia Opossum in an urban environment during the winter.

The Setup

To start, we load the available packages and extract the specific data of interest.

# Load packages
  library(Diel.Niche)
  library(lubridate)
#> 
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:base':
#> 
#>     date, intersect, setdiff, union
  library(coda)
  library(ggplot2)
  library(bayesplot)
#> This is bayesplot version 1.10.0
#> - Online documentation and vignettes at mc-stan.org/bayesplot
#> - bayesplot theme set to bayesplot::theme_default()
#>    * Does _not_ affect other ggplot2 plots
#>    * See ?bayesplot_theme_set for details on theme setting

# Define a year variable
  diel.data$min_year = year(as.POSIXct(diel.data$min_date, format = "%m/%d/%Y"))

# Extract winter data in 2018 for the Virginia Opossum
  winter = subset(diel.data, season=="Winter" & min_year=="2018" & scientificName=="Didelphis virginiana")

# Visualize the data
  head(winter)
#>           scientificName twilight day night trap_nights nsite   min_date
#> 111 Didelphis virginiana        5  20    39        1616   131 12/29/2018
#>      max_date mean_lat mean_lon season       country   phylum    class
#> 111 1/24/2019 41.87236 -87.8423 Winter United States Chordata Mammalia
#>               order      family             Project unit_type      Common_name
#> 111 Didelphimorphia Didelphidae UWIN_Chicago_IL_USA     28day Virginia Opossum
#>     Activity_Literature min_year
#> 111           Nocturnal     2018

For simplicity, we extract the count data and assign it to the object y.

y=data.frame(twilight=winter$twilight,
             day=winter$day, 
             night=winter$night)
rownames(y)=winter$Common_name

y
#>                  twilight day night
#> Virginia Opossum        5  20    39

Model Comparison and Parameter Estimation

We are now ready to compare models using the Traditional hypothesis set, which includes four models: diurnal, nocturnal, crepuscular, and Traditional Cathemeral. We can confirm that this is the hypothesis set of interest by plotting the set together using the plot.diel function as,

plot.diel(hyp=hyp.sets("Traditional"))

To fit our data (y), we use the function ‘diel.fit’, specifying the hypothesis set of interest and MCMC inputs.

out = diel.fit(y = as.matrix(y),hyp.set=hyp.sets("Traditional"),
             post.fit = TRUE, 
             n.chains = 3, n.mcmc=5000,burnin = 1000)
#> Data checks Complete.
#> Calculating Bayes Factors...
#> Posterior Sampling...
#> The most supported model is: 
#>  Cathemeral (Traditional)

The most supported model is Traditional Cathemeral. We can examine the posterior model probabilities for all models as,

out$bf.table
#>    Prior   Posterior
#> D   0.25 0.000000000
#> N   0.25 0.004886813
#> CR  0.25 0.000000000
#> C   0.25 0.995113187

which indicate that we are very certain that Traditional Cathemeral (‘C’) is supported. Notice that without specifying a prior probabilty for each model, it is assumed to be equal among the models.

Checks

As long as we fit our model with multiple chains, we will get an estimate of the Gelman-Rubin convergence diagnostic for our most supported model, as

out$ms.gelm.diag
#> Potential scale reduction factors:
#> 
#>      Point est. Upper C.I.
#> p1_1          1          1
#> p1_2          1          1

Since the point estimate and upper intervals are at one, there is no evidence that the parameters have not converged to their posterior distributions. However, we should also plot the parameter chains to visually inspect this (next section).

We can also examine whether there is evidenve of lack of it betweem our most supported model and the data using a posterior predictive check, which provides a probability of fit (‘ppp’). As long as this value is not near 0 or 1 then there is no evidence of departure of fit to the data.

out$ms.ppc
#>   X2_obs  X2_pred      ppp 
#> 2.006079 1.984158 0.507750

This model indicates that it fits the data with a probability value of NA.

Visualizations

First, let’s plot our parameter chains to confirm that posterior distributions are properly converged.

plot(coda::as.mcmc(out$post.samp.ms.model))

We see that the chains are highly overlapping and thus confirm that posterior distributions have converged.

Next, lets examine the posterior probabilities of activity in the three diel periods,

# Plotting
  plot_title <- ggplot2::ggtitle("Posterior distributions",
                        "with medians and 95% intervals")
  bayesplot::mcmc_areas(out$post.samp.ms.model, prob = 0.95) + plot_title

Alternatively, using base R to plot the same posteriors,

# Plotting
  plot(0,xlim=c(0,1),ylim=c(0,12),ylab="Probability Density", xlab="Probability of Activity",
       main="Virginia Opossum \n Winter 2019")
  polygon(density(out$post.samp.ms.model[,1]), col = "#A73030FF")
  polygon(density(out$post.samp.ms.model[,2]), col = "#EFC000FF")
  polygon(density(out$post.samp.ms.model[,3]), col = "#0073C2FF")
  legend("topright", legend=c("P(twilight)","P(daytime)","P(nighttime)"),lwd=8,
         col=c("#A73030FF","#EFC000FF","#0073C2FF"))

Alternatively, we can use the plot.diel function in the Diel.Niche package that uses plotly to output a 3d plot of the posterior samples on top of the hypothesis set. We execute the ‘plot.diel’ function simply using the generic function ‘plot’ as the objecti out from diel.fit has a unique class identifier (diel).

plot(out)

The black dots are the entire posterior samples plotted within the Traditional Cathemeral hypothesis. This point cloud shows us that the Virginia Opossum is largely active during night and day. To look at the posterior median and 95% credible intervals, we can extract these as,

apply(out$post.samp.ms.model,2,quantile, probs=c(0.025,0.5,0.975))
#>         p_crep_1   p_day_1 p_night_1
#> 2.5%  0.03446237 0.2093689 0.4762338
#> 50%   0.08606358 0.3114301 0.5979032
#> 97.5% 0.16803688 0.4295092 0.7121191

Conclusions

What did we learn? We found that during the winter of 2018 in the urban environment of Chicago, IL that there is strong evidence that the Virginia Opossum is Cathemeral. Specifically, that the Opossum is active most during the nighttime and the daytime.